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ABSTRACT 


For  purposes  of  studying  the  lateral  heterogeneity  as  well  as  for  ultimately  predicting  seismograms  for  this  region, 
we  construct  a  new  3-D  S-velocity  model  by  jointly  inverting  a  variety  of  different  seismic  data.  We  jointly  invert 
regional  waveforms,  surface  wave  group  velocity  measurements,  teleseismic  S  arrival  times,  and  crustal  thickness 
estimates  from  receiver  functions,  refraction  lines,  and  gravity  surveys.  These  data  types  have  complementary 
resolving  power  for  crust  and  mantle  structures,  vertical  and  lateral  variations,  shallow  and  deep  mantle  features, 
local  and  global  structure.  Therefore,  a  joint  inversion  of  these  data  sets  might  help  unravel  the  complexity  of  this 
tectonically  diverse  area.  These  measurements  are  made  from  a  combination  of  mantle  investigation  of  the  deep 
suture  between  Europe  and  Africa  (MIDSEA),  Program  for  Array  Seismic  Studies  of  the  Continental  Lithosphere 
(PASSCAL),  GeoScope,  Geofon,  Global  Seismographic  Network  (GSN),  International  Deployment  of 
Accelerometeres  (IDA),  MedNet,  national  networks,  and  local  deployments  throughout  the  study  region  which 
extends  from  the  western  Mediterranean  region  to  the  Hindu  Kush  and  encompasses  northeastern  Africa,  the 
Arabian  peninsula,  the  Middle  East,  and  part  of  the  Atlantic  Ocean  for  reference.  We  have  fitted  the  waveforms  of 
regional  S  and  Rayleigh  waves  from  over  3800  seismograms  using  Partitioned  Waveform  Inversion.  We  include 
over  3000  crustal  thickness  estimates  from  receiver  functions,  gravity  measurements,  and  refraction  profiles.  We 
include  Rayleigh  wave  group  velocities  for  hundred  thousands  of  paths  transecting  the  region.  We  have  over  3000 
teleseismic  S  arrival  times  measured  through  cross  correlation  and  over  170000  from  picks  originally  reported  to  the 
International  Seismological  Centre  (ISC).  We  demonstrate  our  inversion  methodology  and  discuss  results  from 
combining  these  new  measurements  in  a  joint  inversion  for  a  three-dimensional  S-velocity  model.  Finally,  we  will 
compute  synthetic  seismograms  and  travel-time  surfaces  through  the  new  model. 
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OBJECTIVES 

Our  primary  objective  is  developing  a  new  3-D  ^-velocity  model  for  the  Middle  East  and  Mediterranean  region, 
including  North  Africa,  southern  Europe,  and  Arabia  that 

1)  is  resolved  in  aseismic  regions, 

2)  is  resolved  throughout  the  upper  mantle  (to  660  km), 

3)  resolves  laterally  varying  crustal  thickness, 

4)  contains  laterally  varying  vertical  velocity  gradients, 

5)  is  simultaneously  compatible  with  multiple  data  sets, 

6)  utilizes  several  recent,  unique  waveform  data  sets,  and 

7)  includes  uncertainties  of  the  model  parameters. 

These  features  would  increase  the  model’s  ability  to  predict  and  calibrate  regional  travel  times  and  waveforms, 
thereby  providing  improved  event  locations,  focal  mechanisms,  and  other  event  discriminants. 

Secondly,  we  aim  to  convert  the  3-D  S- velocity  model  to  a  3-D  P- velocity  model,  using  both  literature  on  elastic 
properties  (and  their  partial  derivatives  with  temperature  and  pressure)  of  mantle  rocks  and  empirical  information 
provided  by  measured  arrival  times  of  teleseismic  P  and  .Pms  waves. 

Our  third  objective  is  to  test  both  the  S-and  P-wave  models’  ability  to  predict  regional  P  and  S  travel  times,  deflect 
wave  paths  and  deform  waveforms,  and  assess  their  effects  first  on  the  studied  seismograms  (travel  times  and 
amplitudes)  and  subsequently  on  the  3-D  models  derived  from  these  data. 

This  report  focuses  on  the  primary  objective,  a  new  3-D  S  velocity  model. 
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Figure  1.  Topographic  map  of  the  study  region.  The  pink  line  is  the  NUVEL1-A  (DeMets  et  al.,  1990) 
representation  of  the  Eurasia-Africa-Arabia  plate  boundary.  The  study  region  is  one  of  most 
tectonically  complex  in  the  world  and  is  characterized  by  strong  lateral  variations. 

The  study  region  is  centered  around  the  Africa- Arabia-Eurasia  triple  junction  (Figure  1),  and  extends  west  to  the 
Africa-Eurasia-North  America  junction  at  the  Azores,  just  off  the  map,  and  east  to  the  Arabia-Eurasia-Indian  Plate 
junction.  The  NUVEL1-A  (DeMets  et  al.,  1990)  representation  of  these  plate  boundaries  is  shown  by  the  pink  line  in 
Figure  1.  The  interaction  of  these  six  major  tectonic  plates  with  each  other  and  with  several  microplates  within  an 
area  of  one  quarter  of  the  Earth’s  circumference  yields  this  region  rich  with  tectonic  complexity.  The  three- 
dimensional  structure  of  the  upper  mantle  and  crust  are  correspondingly  complex.  For  example,  while  our  region 
comprises  only  a  few  percent  of  the  Earth’s  surface,  it  is  home  to  more  than  half  of  all  crustal  types  defined  by 
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Mooney  et  al.  (1998)  to  describe  the  Earth’s  crust  globally.  At  the  same  time,  the  region  is  home  to  one  of  the  most 
pervasive,  though  laterally  varying  low-velocity  layers  found  beneath  continental  regions,  which  is  interspersed  with 
high-velocity  bodies,  related  to  Tethyan  subduction,  that  extend  hundreds  of  km  vertically  through  the  upper  mantle 
and  into  the  lower  mantle  (Van  der  Lee,  1990;  Bijwaard  et  al.,  1998;  Marone  et  al.,  2004;  Schmid  et  al.,  2007).  We 
plan  to  capture  various  renditions  of  this  structurally  complicated  part  of  the  world  in  one  S-velocity  model  through 
the  joint  inversion  of  several  different  types  of  seismic  data  simultaneously;  the  new  model  will  refine  our 
understanding  of  the  structure  and  tectonics  in  this  region  of  the  Earth.  The  data  types  we  combine  are  constraints 
from  independent  studies  on  the  depth  to  the  Moho,  fundamental-mode  Rayleigh-wave  group  velocity 
measurements,  waveform  fits  of  regional  S  and  Rayleigh  waves,  and  arrival  times  of  teleseismic  S  waves. 


RESEARCH  ACCOMPLISHED 


Figure  2.  Map  of  raypath  coverage  showing  all  source-receiver  paths  for  the  region  combining  events  (yellow 
circles)  with  Mb  >  5.0  to  all  recording  broadband  stations  (red  triangles)  at  offsets  between  5  and  50 
degrees.  This  represents  approximately  3,800  paths  from  both  MIDSEA  and  LLNL  databases. 

Regional  Waveform  Fits 

We  have  examined  over  13,000  available  waveforms  from  the  Lawrence  Livermore  National  Laboratory  (LLNL) 
database  which  sample  North  Africa  and  the  Middle  East  and  successfully  fit  about  3800  of  them  using  the  non¬ 
linear  inversion  procedure  employed  by  previous  partitioned  waveform  inversion  studies  (Van  der  Lee  and  Nolet, 
1997;  Marone  et  al.,  2004).  The  great-circle  wave  paths  for  these  seismograms  are  shown  in  Figure  2.  We  have 
estimated  path-averaged  S-velocity  structures  for  these  paths  using  the  same  starting  S-velocity  model  but  a  different 
crustal  thickness  (in  5  km  increments)  for  each  path,  based  on  a  priori  reported  estimates  (Van  der  Lee  and  Nolet, 
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1997;  Marone  et  al.,  2004).  The  starting  model  often  predicts  significant  phase  differences  relative  to  the  data  for 
both  the  S-  and  Rayleigh  waves.  The  inversion  procedure  estimates  the  perturbations  to  the  starting  model  by  non¬ 
linear  optimization  (Nolet,  1990;  van  der  Lee  and  Nolet,  1997).  For  the  computation  of  synthetic  waveforms  we  also 
need  to  know  the  source  mechanisms  of  the  events.  We  typically  use  the  centroid  moment  tensors  (CMT)  from  the 
Harvard  or  Global  CMT  project,  and  hypocenters  from  Engdahl’s  reprocessed  ISC  database.  The  frequency  content 
of  the  data  and  synthetic  is  different  for  each  fit  and  generally  falls  within  the  band  0.006  to  0.1  Hz. 

The  optimal  shear- velocity  models  for  some  of  the  paths  in  Figure  2  have  shown:  1)  faster  crustal  velocities  and 
slightly  faster  mantle  velocities  for  the  Nubian  shield,  which  is  consistent  with  work  by  Mokhtar  and  Al-Saeed 
(1994),  Sandvol  et  al.  (1998),  Rodgers  et  al.  (1999),  and  Julia  et  al.  (2003);  2)  lower  velocities  in  the  crust  and  upper 
mantle  for  path  crossing  the  East  African  Rift;  3)  faster  crustal  and  low  sub-Moho  velocities  path  from  the  Gulf  of 
Aden  across  the  Arabian  Shield;  and  4)  faster  velocities  for  paths  traversing  the  Russian  Platform  to  the  north.  It  is 
worth  noting  that  the  Red  Sea  broke  up  the  Nubian- Arabian  Shield  and  these  provinces  could  have  similar  crustal 
petrology.  Both  provinces  have  volcanics  and  it  has  been  speculated  that  mafic  intrusion  may  explain  the  higher 
crustal  velocities.  The  inferred  mantle  velocities  beneath  the  Nubian  Shield  are  slightly  faster  than  the  starting 
model.  However,  reported  mantle  velocities  beneath  the  Arabian  Shield  are  lower  than  average  (e.g.,  Mokhtar  and 
Al-Saeed,  1994;  Rodgers  et  al.,  1999,  Maggi  and  Priestley,  2005). 

The  total  ray  path  coverage  is  extensive  from  both  the  original  MIDSEA  dataset  (Marone  et  al.,  2004)  and  from 
data  in  the  LLNL  database  (shown  in  Figure  2).  We  are  able  to  extend  the  model  area  significantly  to  the  east  and 
achieve  dense  sampling  of  the  Arabian  Shield,  Iran,  Afghanistan,  and  Pakistan.  We  also  have  good  coverage  of 
northeast  Africa,  the  Red  Sea,  and  parts  of  the  Russian  Platform.  We  combine  the  1-D,  path-averaged,  constraints 
obtained  from  these  -3,800  waveform  fits  with  other  data  sets  in  the  joint  inversion  for  shear  velocity  structure. 


Surface  Wave  Group  Velocities 


Figure  3.  Total  wave  path  coverage  for  the  Rayleigh-wave  group  velocity  data  set.  Yellow  dots  are 

earthquakes,  and  red  triangles  are  seismic  stations.  White  lines  are  great-circle  wave  paths. 


We  have  measured  group  velocities  of  fundamental-mode  Rayleigh  waves  recorded  at  MID  SEA  and  other  stations 
in  the  region  and  used  them  to  update  previous  group  velocity  maps  (e.g.,  Pasyanos,  2005).  Figure  3  shows  the 
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great-circle  wave  paths  for  which  we  included  fundamental-mode  Rayleigh  wave  group  velocities  in  our  joint 
inversion.  Compared  to  other  data  types  we  include,  lateral  coverage  is  best  for  this  group  velocity  data  set,  though 
vertical  coverage  is  provided  mainly  by  the  teleseismic  arrival  times  (Figure  5)  and  regional  S  and  Rayleigh 
waveforms  (Figure  2). 

Pasyanos  (2005),  for  example,  used  these  group  velocities  in  a  typical  two-step  inversion  process,  first  inverting  for 
lateral  variations  in  group  velocity  for  a  given  frequency  of  the  Raleigh  wave  and  then  inverting  the  group  velocities 
for  different  frequencies  at  a  given  locality  for  a  local  shear- velocity  model.  To  allow  a  self-consistent  treatment  of 
all  data  in  our  joint  inversion  we  modified  these  steps  such  that  first  we  express  the  measured  group  velocities  as 
path- averages  and  integrals  over  the  S-velocity  distribution  with  depth.  For  each  path  we  thus  obtain  as  many 
constraints  on  the  laterally  averaged  S-velocity  as  we  have  measurements  at  different  periods.  Typically  these 
equations  have  a  high  degree  of  redundancy  with  respect  to  describing  the  S-velocity  distribution  with  depth  and  are 
far  from  independent.  We  distill  the  most  powerful  constraints  out  of  each  such  set  of  equations  using  a  singular 
value  decomposition  to  recombine  the  equations  into  independent  ones  and  discarding  constraints  associated  with 
small  singular  values.  This  process  also  provides  meaningful  error  estimates  that  we  later  use  to  weight  the 
constraints  with  smaller  error  estimates  more  strongly  during  the  joint  inversion.  We  then  solve  these  subsets  of 
equations  for  the  3-D  S-velocity  distribution  by  combining  all  of  them,  with  or  without  constraints  from  our  other 
data  sets,  in  a  single,  linear,  least-squares,  joint  inversion.  We  benchmarked  this  new  2-step  procedure  against 
S-velocity  models  obtained  with  the  traditional  2-step  procedure  and  find  that  it  performs  well  and  efficiently.  The 
added  advantage  is  that  we  can  jointly  and  self-consistently  invert  the  group  velocity  data  together  with  our  other 
data. 

Crustal  Thickness  Constraints 
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Figure  4.  Total  point  coverage  for  independent  constraints  on  Moho  depth  used  in  our  joint  inversion.  The 
points  are  color-coded  with  respect  to  a  priori  determined  Moho  depths. 

While  the  fundamental-mode  Rayleigh-wave  group  velocities  and  regional  S  and  Rayleigh  waveforms  have 
significant  sensitivity  to  Moho  depth,  they  cannot  uniquely  resolve  Moho  depth.  We  therefore  include  independent 
estimates  of  crustal  thickness  as  point  constraints  in  the  joint  inversion  and  thus  compile  such  measurements  from  a 
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large  number  of  published  studies.  A  partial  list  of  these  studies  is  provided  in  Marone  et  al.  (2003),  and  part  of 
these  points  stem  from  the  database  used  for  CRUST5.1.  We  have  added  new  estimates  of  Moho  depth  from  more 
recent  studies  and  interactively  resolved  or  removed  conflicting  data  and  outliers  from  the  data  set.  The  majority  of 
this  data  set  is  from  receiver  function  studies.  Remaining  points  are  from  active-source  studies  and  some  are  from 
gravity  surveys.  The  active-source  constraints  were  assigned  the  lowest  error  and  the  Moho-depth  estimates  from 
gravity  were  assigned  the  largest  errors.  For  the  oceans  we  use  a  constraint  of  10  km  for  Moho  depth,  but  only  for 
points  also  covered  by  data  from  our  other  data  sets. 


Teleseismic  N-wave  Arrival  Times 


Figure  5.  Earthquakes  (circles)  and  stations  (triangles)  for  which  we  have  teleseismic  S-wave  travel-times  in 
our  joint  inversion.  Blue  stations  and  cyan  events  are  used  for  cross-correlation  measurements  and 
red  stations  and  yellow  events  are  from  Engdahl’s  (Engdahl  et  al.,  1998)  database  of  arrival  picks. 

Figure  5  shows  a  map  of  ah  stations  and  earthquakes  that  produced  S-wave  arrival  times  that  we  use  in  our  joint 
inversion.  The  majority  of  these  are  from  Engdahl’ s  reprocessed  ISC  data  set  spanning  38  years  of  recording. 
Nevertheless,  because  of  their  heterogeneity,  these  arrival  times  were  assigned  a  higher  error  than  those  we 
measured  through  cross-correlation.  The  waveforms  of  teleseismic  S  recorded  at  MIDSEA  and  other  stations  in  the 
study  region  were  cross-correlated  with  the  method  of  VanDecar  and  Crosson  (1990),  yielding  high-accuracy 
relative  S-wave  arrival  times.  Figure  6  shows  those  portions  of  the  teleseismic  S-wave  paths  in  our  data  set  that  he 
within  a  half  degree  from  the  cross  section.  Both  cross  sections  show  large  data  gaps  in  the  upper  mantle  beneath 
Arabia  and  North  Africa,  illustrating  the  rationale  for  jointly  inverting  these  data  with  surface  wave  data,  which  are 
sensitive  to  upper-mantle  structure  ah  along  their  way  from  event  to  station.  The  cross  sections  in  Figure  6  show  that 
the  top  of  the  lower  mantle  is  well  covered  by  these  data  and  that  the  ray  paths  are  near- vertical  in  the  uppermost 
mantle,  providing  excellent  lateral  resolving  power.  This  lateral  resolving  power  excellently  complements  the 
vertical  resolving  power  of  the  regional  S  and  Rayleigh  waveform  fits.  Resolving  Moho  depth  from  uppermost 
mantle  and  crust  velocities  is  aided  by  the  independent  constraints  on  the  Moho  and  the  fundamental-mode  group 
velocity  measurements. 
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Figure  6.  A  NS  and  EW  cross  section  through  the  study  region  showing  segments  of  ray  paths  of  teleseismic  S 
waves  used  in  our  joint  inversion  that  lie  within  half  a  degree  from  the  cross  section. 


Inversion  Results 
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Figure  7.  Moho  map  from  joint  inversion  of  Moho  constraints  (Figure  4),  regional  waveform  fits  (Figure  2), 
Rayleigh  wave  group  velocities  (Figure  3),  and  teleseismic  arrival  times  (Figure  5). 

The  benefit  of  our  joint  inversion  lies  in  the  fact  that  the  different  types  of  data  sets  we  use  cover  the  same  general 
volume  of  the  mantle  and  crust  along  the  Tethys  suture  but  cover  the  volume  in  very  different,  but  complementary 
ways;  our  data  sets  are  both  redundant  and  complementary.  The  redundancy  is  needed  to  increase  accuracy  and  to 
ensure  that  all  data  sets  measure  the  same  structural  phenomena.  The  data  sets  need  to  be  complementary  to  increase 
resolving  power  over  a  larger  volume  of  mantle  and  crust  and  thereby  reduce  trade-offs  inherent  in  each  type  of 
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seismic  data  set.  Combining  the  data  sets  therefore  not  only  enhances  resolving  power  where  the  region  is  covered 
by  more  than  one  data  sets,  but  resolution  is  also  enhanced  in  the  regions  covered  by  only  one  data  set  (such  as  the 
lower  mantle  or  the  upper  crust)  because  of  the  reduced  trade-offs  with  other  regions  of  the  model. 

We  have  developed  software  that  handles  the  joint  inversion  of  constraints  from  regional  waveform  fits,  crustal 
estimates,  group  velocities,  and  teleseismic  arrival  times.  An  early  version  of  the  joint  inversion  code  has  been  tested 
on  the  teleseismic  S  arrival  time  data  set  of  Schmid  et  al.  (2006)  and  the  data  derived  from  regional  waveform  fitting 
from  Marone  et  al.  (2004).  The  results  show  only  a  percent  or  two  increase  in  the  variance  reduction  obtained  in  the 
linear  inversion  of  both  data  sets  compared  to  their  individual  inversions  (Schmid  et  al.,  2006).  The  resolving  power 
of  the  joint  data  sets,  however,  has  increased  dramatically:  the  teleseismic  data  add  more  lateral  resolution  to  the 
regional  waveform  data,  while  the  regional  waveform  data  add  more  depth  resolution  to  the  teleseismic  data.  The 
resolved  depth  range  of  the  combined  data  sets  has  doubled  with  respect  to  their  individual  depth  ranges.  The 
regional  waveform  data  resolve  the  upper  mantle  more  strongly  while  the  teleseismic  arrival  times  resolve  the  lower 
mantle  more  strongly.  The  resolving  power  of  the  combined  data  is  superior  to  that  of  each  of  the  data  sets  alone. 

The  joint  inversion  reduces  the  variance  in  the  data  sets  by  40  %  for  the  teleseismic  delay  times,  50  %  for  the 
Rayleigh-wave  group  velocities,  75  %  for  the  Moho  point  constraints,  and  90  %  for  the  regional  waveform  fits. 

Figure  7  shows  the  joint  inversion  result  for  Moho  depth.  The  map  shows  a  good  resemblance  with  the  point 
constraints  in  Figure  4,  but  is  much  smoother  because  the  regional  waveforms  and  group  velocity  data,  as  well  as 
regularization  constraints  in  the  joint  inversion  provide  a  smooth  interpolation  between  the  point  constraints.  The 
map  is  also  broadly  consistent  with  CRUST2.0,  except  in  northern  Africa,  where  the  crust  from  our  joint  inversion  is 
about  5  km  thinner.  Resolution  tests  indicate  that  our  combined  data  sets  can  resolve  the  crustal  thicknesses  of 
CRUST2.0. 


-300-233-166-100  -33  33  100  166  233  300  m/s 

Figure  8.  S-velocity  perturbations  relative  to  1-D  model  MEAN  (Marone  et  al.,  2004)  at  a  depth  of  100  km 
(left)  and  700  km  (right)  from  our  joint  inversion. 


The  left  panel  of  Figure  8  shows  upper  mantle  shear  velocity  anomalies  from  the  joint  inversion  at  a  depth  of 
100  km.  The  Mid- Atlantic  Ridge,  Algero-Provencal  and  Pannonian  Basins,  Tyrrhenian  and  Aegean  Seas,  Arabian 
and  Nubian  shields,  and  the  Anatolian  and  Iranian  Plateaus  are  underlain  by  significant  low  velocities.  High 
velocities  are  imaged  beneath  the  West- African  Craton,  to  the  northeast  of  the  Tomquist-Tesseire  Zone,  beneath  the 
Hellenic  and  Calabrian  Arcs,  and  beneath  the  Arabian  Platform  and  Zagros  Fold  Belt.  The  right  panel  of  Figure  8 
shows  that  the  anomalies  get  weaker  below  the  upper  mantle  and  are  dominated  by  high-velocity  anomalies  in  the 
top  of  the  lower  mantle.  These  high-velocity  anomalies  appear  to  form  an  EW  belt,  roughly  parallel  to  the  Tethys 
suture. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


The  broad  consistency  between  seismic  velocity  anomalies  inferred  from  teleseismic  arrival  times,  Rayleigh  wave 
group  velocities,  and  regional  waveforms  shown  here  implies  that  these  different  types  of  data  sets  have  recorded 
the  same  structural  phenomena,  despite  differences  in  size  and  character  between  typical  sensitivity  kernels  for  each 
data  set.  This  conclusion  is  further  supported  by  an  analysis  of  how  teleseismic  delay  times  depend  on  frequency 
(Schmid  et  ah,  2006)  and  that  the  teleseismic  arrival  times  and  the  regional  waveforms  are  highly  complementary. 
The  shared  sensitivity,  though  different  in  character,  of  receiver  functions  and  Rayleigh  wave  group  velocities  to 
crustal  structure  helps  to  separate  crustal  effects  on  the  observed  data  from  mantle  effects,  when  included  in  the  joint 
inversion. 

Through  a  joint  inversion  of  teleseismic  S-wave  arrival  time  delays,  waveform  fits  of  regional  S-  and  Rayleigh 
waves,  group  velocity  measurements  of  fundamental-mode  Rayleigh  waves,  and  independent  constraints  on  Moho 
depth,  we  have  achieved  considerable  variance  reduction  in  each  data  set  simultaneously.  The  3-D  S -velocity  model 
resulting  from  this  joint  inversion  shows  high  and  low- velocity  anomalies  that  are  broadly  consistent  with  previous 
studies,  restricted  to  a  smaller  region  or  a  more  limited  data  set.  Our  new  model  may  be  converted  to  a  P-velocity 
model  and  used  for  locating  events. 

Our  final  steps  in  this  project  are  to  carry  out  this  conversion  to  P-velocities  and  computing  arrival  time  surfaces  and 
synthetic  seismograms  for  the  purpose  of  potentially  improving  event  location  procedures. 
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